function [alpha,beta,perc]=see_density(mu,sd,density,lb,ub,varargin)
%Function creates plot of density along with plot of the log of density.
%   Densities Supported: Beta,Gamma,Normal,Inverse Gamma 1, Weibull
%   Inputs: 
%     mu: mean of distribution
%     sd: standard dev of distribution
%     lb,ub: lower and upper bound for plotting
%     varargin: An option for the number of draws to use when computing the
%               percentiles of the IG1 distribution
%   Outputs: 
%       alpha,beta: The natural parameters of the distribution
%       perc: The percentiles of the distribution
%   Notes: In the case of IG1 [alpha,beta]=[v,s]
%          In the case of weibull [alpha,beta]=[lambda,k]
%          In the case of Normal, [alpha,beta]=[mu,sd]

    step=.001; %gridsize for plotting the distribution
    switch density
        
        case {'beta','Beta'}
            %Find the natural parameters
            if mu > 1; error('Mean must be less than 1'); end
            
            beta=mu*( (1-mu)^2 )/(sd^2) - ( 1 - mu );
            alpha=mu*beta/(1-mu); 
            
            if beta < 0 || alpha < 0 
                error('Check the settings of your beta') 
            end
            
            %% Plot the density
            subplot(2,1,1)
            [percentiles junk] = plot_beta(alpha,beta,lb,step,ub);
            
            axis tight; 
            
            %Grab the values from the plot, which will be needed in
            %computing the second plot
            data=findall(gcf,'type','line');
            xData=get(data,'xdata');
            yData=get(data,'ydata');
            legend off
            
            %Plotting the percentiles: Note, in order to get the
            %percentiles to appear in a legend, we need to put them on the
            %plot.  I chose white so the points wouldn't show up
            perc = cell2num(percentiles);
            perc=perc(3:end,:);
            hold on
            for k=1:size(perc,1)
                plot(xData{3}(1),yData{3}(1),'Color','white')
                legendLabels{k}=[num2str(perc(k,1)),' = ', num2str(perc(k,2))];
                hold on
            end
            legHand=legend(legendLabels,'location','Best','Orientation','vertical');
            
            %Plot the log of the density
            subplot(2,1,2)
            plot(xData{3},log(yData{3}))
            title('Log Density of Beta')
            
        case {'gamma','Gamma'}
            beta=(sd^2)/mu;
            alpha=mu/beta;

            [percentiles,junk] = plot_gamma(alpha,beta,lb,step,ub);
            data=findall(gcf,'type','line');
            xData=get(data,'xdata');
            yData=get(data,'ydata');
            close
            subplot(2,1,1)
            plot(xData{3},yData{3})
            titleString=['Gamma(',num2str(alpha),',',num2str(beta),')  ','mean: ',num2str(mu),' s.d.: ', num2str(sd)];
            title(titleString)
            perc = cell2num(percentiles);
            perc=perc(3:end,:);
            hold on
            for k=1:size(perc,1)
                plot(lb,lb,'Color','white')
                legendLabels{k}=[num2str(perc(k,1)),' = ', num2str(perc(k,2))];
                hold on
            end
            legHand=legend(legendLabels,'location','Best','Orientation','vertical');

            subplot(2,1,2)
            plot(xData{3},log(yData{3}))
            title('Log Density of Gamma')   
            
        case {'normal', 'Normal'}
            
            subplot(2,1,1)
            grid=(lb:step:ub);
            plot(grid,normpdf(grid,mu,sd))
            titleString=['Normal(',num2str(mu),',',num2str(sd),')  ','mean: ',num2str(mu),' s.d.: ', num2str(sd)];
            title(titleString)
          
            percentiles=[.001,.01,.05,.1,.5,.9,.95,.99,.999]';
            perc=norminv(percentiles,mu,sd);
            hold on
            for k=1:size(perc,1)
                plot(mu,0,'Color','white')
                legendLabels{k}=[num2str(percentiles(k)),' = ', num2str(perc(k))];
                hold on
            end
            legHand=legend(legendLabels,'location','Best','Orientation','vertical');

            subplot(2,1,2)
            plot(grid,log(normpdf(grid,mu,sd)))
            title('Log Density of Normal')
            alpha=mu;
            beta=sd;
        
        case {'IG1','ig1','InverseGamma1','IGone','igone','Inverse Gamma 1'}
        
            if nargin>5
                ndr=varargin{1};
            else
                ndr=1E5;
            end
            
            pvec=[0.001 0.01 0.05 0.1 0.5 0.9 0.95 0.99 0.995]; 
            if sd > 40 
                sd=inf;
                display('SD set to infinity')
            end 
            [s,v]=inverse_gamma_specification( mu , sd ,1 ); 
            gr=linspace(lb,ub);
            pd=zeros(100,1); 
            for ii=1:100; 
                pd(ii)=pdf_igone(gr(ii),v,s) ; 
            end 
            mod=( (v/(v+1))^0.5 )*sqrt(s/v); 
    
            dr=zeros(ndr,1); 
            sim=1./gamrnd(v/2,2/s,ndr,1); 
            sim=sqrt(sim);
            pvec=floor(ndr*pvec); 
            sims=sort(sim);
            perc=[(pvec/ndr)' sims(pvec)]; 
            
            subplot(2,1,1)
            plot(gr,pd)
            titleString=['PDF IG1 Mean ',num2str(mu),' SD: ',num2str(sd), ' ; v ', num2str(v) , ' s= ',num2str(s)];
            title(titleString)
            hold on
            for k=1:size(perc,1)
                plot(gr(1),pd(1),'Color','white')
                legendLabels{k}=[num2str(perc(k,1)),' = ', num2str(perc(k,2))];
                hold on
            end
            
            legHand=legend(legendLabels,'location','Best','Orientation','vertical');
            
            subplot(2,1,2)
            plot(gr,log(pd))
            title('Log Density of Inverse Gamma 1')
            alpha=v;
            beta=s;
        
        case {'weibull','Weibull'}
            display('Solving for natural parameters')
                   
            weibFunc = @(k) gamma(1+1./k).^2./gamma(1+2./k)-mu^2/(sd^2+mu^2);
            grid=(1E-3:.1:100);
            % Find an approximate solution
            [val loc] = min(abs(feval(weibFunc,grid)));
            % Find a more exact solution
            [beta fval exit_flag]=fzero(weibFunc,grid(loc));
        
            if exit_flag ~= 1
                error('Failed to solve for weibull parameters');
            end

            alpha = mu/gamma(1+1/beta);
            
            if alpha <= 0 ||  beta < 0
                error('The mean/SD values lead to invalid parameters')
            end
                
            pvec=[0.001 0.01 0.05 0.1 0.5 0.9 0.95 0.99 0.995]';
            perc=wblinv(pvec,alpha,beta);

            subplot(2,1,1)
            gr=linspace(lb,ub);
            pd=wblpdf(gr,alpha,beta);
            plot(gr,pd)
            titleString=['PDF Weibull(',num2str(alpha),',',num2str(beta),')   mean: ', num2str(mu), '  sd: ', num2str(sd)];
            title(titleString)
            hold on
            for k=1:length(perc)
                plot(gr(1),pd(1),'Color','white')
                legendLabels{k}=[num2str(pvec(k)),' = ', num2str(perc(k))];
                hold on
            end
            
            legHand=legend(legendLabels,'location','Best','Orientation','vertical');
            subplot(2,1,2)
            
            plot(gr,log(pd))
            title('Log Density of Weibull')
            perc=[pvec perc];
        
       otherwise
            display('Not a supported density: Try "beta","gama","normal","igone",or "weibull" ')
            
    end
        



end